Fit spatiotemporal diet models with covariates
Max Lindmark
25 November, 2022
Intro
In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.
Load packages & source functions
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(viridis)
library(sdmTMB)
library(mapplots)
library(sp)
library(raster)
library(ggeffects)
library(modelr)
# Source code for lan_lot to UTM
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/lon_lat_utm.R")
# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")
# Trim the map plot to better fit stomach data
xmin2 <- 303000
xmax2 <- 916000
xrange <- xmax2 - xmin2
ymin2 <- 5980000
ymax2 <- 6580000
yrange <- ymax2 - ymin2
theme_set(theme_plot())
plot_map_fc2 <- plot_map_fc +
theme(legend.position = "bottom") +
xlim(xmin2*1.5, xmax2*0.9) +
ylim(ymin2*1.025, ymax2*0.98268)
# Continuous colors
options(ggplot2.continuous.colour = "viridis")Read stomach data
small_cod_stomach <- read_csv("data/clean/small_cod_stomach.csv") %>%
drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>%
rename(density_cod_small = scod_kg_km2,
density_cod_large = lcod_kg_km2,
density_flounder = fle_kg_km2) %>%
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
ices_rect = as.factor(ices_rect),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
temp_sc = (temp - mean(temp)) / sd(temp),
depth_sc = (depth - mean(depth)) / sd(depth),
density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
haul_id = as.factor(haul_id),
species2 = "Small cod")
#> Rows: 1153 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (24): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#>
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#>
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#>
#> converted 'haul_id' from character to factor (0 new NA)
#>
#> new variable 'fyear' (factor) with 6 unique values and 0% NA
#>
#> new variable 'fquarter' (factor) with 2 unique values and 0% NA
#>
#> new variable 'oxy_sc' (double) with 152 unique values and 0% NA
#>
#> new variable 'temp_sc' (double) with 152 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 64 unique values and 0% NA
#>
#> new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#>
#> new variable 'density_cod_large_sc' (double) with 157 unique values and 0% NA
#>
#> new variable 'density_flounder_sc' (double) with 159 unique values and 0% NA
#>
#> new variable 'species2' (character) with one unique value and 0% NA
large_cod_stomach <- read_csv("data/clean/large_cod_stomach.csv") %>%
drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>%
rename(density_cod_small = scod_kg_km2,
density_cod_large = lcod_kg_km2,
density_flounder = fle_kg_km2) %>%
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
ices_rect = as.factor(ices_rect),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
temp_sc = (temp - mean(temp)) / sd(temp),
depth_sc = (depth - mean(depth)) / sd(depth),
density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
haul_id = as.factor(haul_id),
species2 = "Large cod")
#> Rows: 2153 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (24): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#>
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#>
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#>
#> converted 'haul_id' from character to factor (0 new NA)
#>
#> new variable 'fyear' (factor) with 6 unique values and 0% NA
#>
#> new variable 'fquarter' (factor) with 2 unique values and 0% NA
#>
#> new variable 'oxy_sc' (double) with 164 unique values and 0% NA
#>
#> new variable 'temp_sc' (double) with 164 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 67 unique values and 0% NA
#>
#> new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#>
#> new variable 'density_cod_large_sc' (double) with 171 unique values and 0% NA
#>
#> new variable 'density_flounder_sc' (double) with 171 unique values and 0% NA
#>
#> new variable 'species2' (character) with one unique value and 0% NA
flounder_stomach <- read_csv("data/clean/flounder_stomach.csv") %>%
drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>%
rename(density_cod_small = scod_kg_km2,
density_cod_large = lcod_kg_km2,
density_flounder = fle_kg_km2) %>%
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
ices_rect = as.factor(ices_rect),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
temp_sc = (temp - mean(temp)) / sd(temp),
depth_sc = (depth - mean(depth)) / sd(depth),
density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
haul_id = as.factor(haul_id),
species2 = "Flounder")
#> Rows: 2579 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (24): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#>
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#>
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#>
#> converted 'haul_id' from character to factor (0 new NA)
#>
#> new variable 'fyear' (factor) with 6 unique values and 0% NA
#>
#> new variable 'fquarter' (factor) with 2 unique values and 0% NA
#>
#> new variable 'oxy_sc' (double) with 153 unique values and 0% NA
#>
#> new variable 'temp_sc' (double) with 153 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 64 unique values and 0% NA
#>
#> new variable 'density_cod_small_sc' (double) with 148 unique values and 0% NA
#>
#> new variable 'density_cod_large_sc' (double) with 160 unique values and 0% NA
#>
#> new variable 'density_flounder_sc' (double) with 161 unique values and 0% NA
#>
#> new variable 'species2' (character) with one unique value and 0% NA
# Load cache
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/main_diet_analysis_cache/html")Read the prediction grid (scales variables wrt data!)
# I need predicted cod and flounder densities for this
pred_grid_1 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(1_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid_2 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(2_2).csv")
#> Rows: 296520 Columns: 18
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): substrate, ices_rect
#> dbl (16): X, Y, depth, year, lon, lat, quarter, oxy, temp, sal, sub_div, den...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- bind_rows(pred_grid_1, pred_grid_2)
# # Flounder
# pred_grid_fle <- pred_grid %>%
# filter(year %in% unique(flounder_stomach$year) & ices_rect %in% unique(flounder_stomach$ices_stat_rec)) %>%
# # scale pred_grid wrt data
# mutate(fyear = as.factor(year),
# temp_sc = (temp - mean(flounder_stomach$temp)) / sd(flounder_stomach$temp),
# oxy_sc = (oxy - mean(flounder_stomach$oxy)) / sd(flounder_stomach$oxy),
# depth_sc = (depth - mean(flounder_stomach$depth)) / sd(flounder_stomach$depth),
# density_cod_large_sc = (density_cod_large - mean(flounder_stomach$density_cod_large)) / sd(flounder_stomach$density_cod_large),
# density_cod_small_sc = (density_cod_small - mean(flounder_stomach$density_cod_small)) / sd(flounder_stomach$density_cod_small),
# density_flounder_sc = (density_fle - mean(flounder_stomach$density_flounder)) / sd(flounder_stomach$density_flounder))
#
# plot_map_fc +
# geom_point(data = pred_grid_fle, aes(X, Y)) +
# facet_wrap(~year)
#
# # Cod
# pred_grid_large_cod <- pred_grid %>%
# filter(year %in% unique(large_cod_stomach$year) & ices_rect %in% unique(large_cod_stomach$ices_stat_rec)) %>%
# mutate(fyear = as.factor(year),
# temp_sc = (temp - mean(large_cod_stomach$temp)) / sd(large_cod_stomach$temp),
# oxy_sc = (oxy - mean(large_cod_stomach$oxy)) / sd(large_cod_stomach$oxy),
# depth_sc = (depth - mean(large_cod_stomach$depth)) / sd(large_cod_stomach$depth),
# density_cod_large_sc = (density_cod_large - mean(large_cod_stomach$density_cod_large)) / sd(large_cod_stomach$density_cod_large),
# density_cod_small_sc = (density_cod_small - mean(large_cod_stomach$density_cod_small)) / sd(large_cod_stomach$density_cod_small),
# density_flounder_sc = (density_fle - mean(large_cod_stomach$density_flounder)) / sd(large_cod_stomach$density_flounder))
#
# plot_map_fc +
# geom_point(data = pred_grid_large_cod, aes(X, Y)) +
# facet_wrap(~year)
#
#
# pred_grid_small_cod <- pred_grid %>%
# filter(year %in% unique(small_cod_stomach$year) & ices_rect %in% unique(small_cod_stomach$ices_stat_rec)) %>%
# mutate(fyear = as.factor(year),
# temp_sc = (temp - mean(small_cod_stomach$temp)) / sd(small_cod_stomach$temp),
# oxy_sc = (oxy - mean(small_cod_stomach$oxy)) / sd(small_cod_stomach$oxy),
# depth_sc = (depth - mean(small_cod_stomach$depth)) / sd(small_cod_stomach$depth),
# density_cod_large_sc = (density_cod_large - mean(small_cod_stomach$density_cod_large)) / sd(small_cod_stomach$density_cod_large),
# density_cod_small_sc = (density_cod_small - mean(small_cod_stomach$density_cod_small)) / sd(small_cod_stomach$density_cod_small),
# density_flounder_sc = (density_fle - mean(small_cod_stomach$density_flounder)) / sd(small_cod_stomach$density_flounder))
#
# plot_map_fc +
# geom_point(data = pred_grid_small_cod, aes(X, Y)) +
# facet_wrap(~year)Plot data in space
full_dat <- bind_rows(small_cod_stomach, large_cod_stomach, flounder_stomach) %>%
pivot_longer(c(tot_feeding_ratio, saduria_feeding_ratio, benthic_feeding_ratio),
names_to = "prey_group",
values_to = "feeding_ratio")
#> pivot_longer: reorganized (tot_feeding_ratio, benthic_feeding_ratio, saduria_feeding_ratio) into (prey_group, feeding_ratio) [was 5885x41, now 17655x40]
## Total feeding ratio
# Q1
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 1 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Total feeding ratio", subtitle = "Quarter 1") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining
# Q4
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 4 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Total feeding ratio", subtitle = "Quarter 4") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining
## Benthic feeding ratio
# Q1
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 1 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Benthic feeding ratio", subtitle = "Quarter 1") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining
# Q4
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 4 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Benthic feeding ratio", subtitle = "Quarter 4") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining
## Saduria feeding ratio
# Q1
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 1 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Saduria feeding ratio", subtitle = "Quarter 1") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining
# Q4
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 4 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Saduria feeding ratio", subtitle = "Quarter 4") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining
Plot data (response vs explanatory variables)
ggplot(full_dat, aes(depth, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(oxy, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(temp, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(density_flounder_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(density_cod_large_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(density_cod_small_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Fit spatiotemporal model to feeding ratio
sdmTMB models with densities as covariates to stomach content data. First make mesh:
# Small cod
small_cod_spde <- make_mesh(small_cod_stomach,
c("X", "Y"),
cutoff = 15,
type = "kmeans",
seed = 42)
plot(small_cod_spde)
# Large cod
large_cod_spde <- make_mesh(large_cod_stomach,
c("X", "Y"),
n_knots = 75,
type = "kmeans",
seed = 42)
plot(large_cod_spde)
# Flounder
flounder_spde <- make_mesh(flounder_stomach,
c("X", "Y"),
n_knots = 75,
type = "kmeans",
seed = 42)
plot(flounder_spde)
# s_cod_tot0 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# tidy(s_cod_tot0)
# sanity(s_cod_tot0)
#
# # Plot residuals in space
# small_cod_stomach$tot_resid0 <- residuals(s_cod_tot0)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid0), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
#
# s_cod_tot <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# tidy(s_cod_tot)
# sanity(s_cod_tot)
#
# # Plot residuals in space
# small_cod_stomach$tot_resid <- residuals(s_cod_tot)
#
# # plot_map_fc2 +
# # geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# # shape = 21, color = "grey50", stroke = 0.2) +
# # facet_wrap(~ year) +
# # scale_fill_gradient2()
# #
# # plot_map_fc2 +
# # geom_point(data = filter(small_cod_stomach, tot_resid < 5), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# # shape = 21, color = "grey50", stroke = 0.2) +
# # facet_wrap(~ year) +
# # scale_fill_gradient2()
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
# # Fit model now with spatial random effect
# s_cod_tot2 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "on",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot2)
# tidy(s_cod_tot2)
# sanity(s_cod_tot2)
#
# small_cod_stomach$tot_resid2 <- residuals(s_cod_tot2)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid2), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
# # Fit model now with spatiotemporal random effect
# s_cod_tot3 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "IID",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot3)
# tidy(s_cod_tot3)
# sanity(s_cod_tot3)
#
# small_cod_stomach$tot_resid3 <- residuals(s_cod_tot3)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid3), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
# # Fit model now with spatiotemporal AND spatial random effect
# s_cod_tot4 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "IID",
# spatial = "on",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot4)
# tidy(s_cod_tot4)
# sanity(s_cod_tot4)
#
# small_cod_stomach$tot_resid4 <- residuals(s_cod_tot4)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid4), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
#
# # Just random effects...
# small_cod_stomach$ices_rect <- as.factor(small_cod_stomach$ices_rect)
#
# s_cod_tot5 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|ices_rect) + (1|haul_id),
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot5)
# tidy(s_cod_tot5)
# sanity(s_cod_tot5)
#
# small_cod_stomach$tot_resid5 <- residuals(s_cod_tot5)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid5), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
#
# # Plot together
# long <- small_cod_stomach %>%
# pivot_longer(c(tot_resid0, tot_resid, tot_resid2, tot_resid3, tot_resid4, tot_resid5), names_to = "Model", values_to = "Residuals")
#
# plot_map_fc2 +
# geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
# shape = 21, color = "grey50", stroke = 0.1) +
# scale_fill_gradient2() +
# facet_wrap(~Model, ncol = 3)
#
# plot_map_fc2 +
# geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
# shape = 21, color = "grey50", stroke = 0.1) +
# scale_fill_gradient2() +
# facet_grid(year ~ Model)
#
#
# # Plot spatial and spatiotemporal random effect
# pred_grid_small <- pred_grid %>%
# filter(quarter == 4) %>%
# mutate(X = X/1000,
# Y = Y/1000) %>%
# filter(year %in% unique(small_cod_stomach$year)) %>%
# filter(X < max(small_cod_stomach$X)) %>%
# filter(X > min(small_cod_stomach$X)) %>%
# filter(Y < max(small_cod_stomach$Y)) %>%
# filter(Y > min(small_cod_stomach$Y))
#
# pred_grid_small <- pred_grid_small %>%
# mutate(fyear = factor(year),
# fquarter = factor(quarter),
# temp_sc = 0,
# depth_sc = 0,
# density_cod_small_sc = 0,
# density_cod_large_sc = 0,
# density_flounder_sc = 0,
# oxy_sc = 0)
#
# spat_pred <- predict(s_cod_tot2, newdata = pred_grid_small)
#
# plot_map_fc2 +
# geom_raster(data = spat_pred, aes(X*1000, Y*1000, fill = omega_s)) +
# scale_fill_gradient2()
#
# spat_pred2 <- predict(s_cod_tot3, newdata = pred_grid_small)
#
# plot_map_fc2 +
# geom_raster(data = spat_pred2, aes(X*1000, Y*1000, fill = epsilon_st)) +
# scale_fill_gradient2() +
# facet_wrap(~year)
#
# AIC(s_cod_tot0, s_cod_tot, s_cod_tot2, s_cod_tot3, s_cod_tot4, s_cod_tot5)Small cod
# Total
s_cod_tot <- sdmTMB(
data = small_cod_stomach,
formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(s_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc + density_cod_large_sc + density_flounder_sc *
#> Formula: oxy_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -5.27 0.41
#> fyear2016 -5.42 0.26
#> fyear2017 -4.72 0.20
#> fyear2018 -5.57 0.21
#> fyear2019 -4.87 0.37
#> fyear2020 -5.59 0.19
#> fquarter4 -0.09 0.31
#> temp_sc -0.18 0.11
#> depth_sc -0.31 0.14
#> density_cod_small_sc -0.11 0.23
#> density_cod_large_sc 0.17 0.20
#> density_flounder_sc 0.04 0.10
#> oxy_sc -0.12 0.16
#> density_flounder_sc:oxy_sc -0.22 0.12
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.77
#>
#> Dispersion parameter: 0.18
#> Tweedie p: 1.59
#> ML criterion at convergence: -2952.284
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_tot)
#> term estimate std.error
#> 1 fyear2015 -5.27421186 0.41041126
#> 2 fyear2016 -5.42038244 0.26174994
#> 3 fyear2017 -4.71694581 0.19506241
#> 4 fyear2018 -5.56547652 0.20792756
#> 5 fyear2019 -4.87241401 0.37423092
#> 6 fyear2020 -5.58949458 0.18710427
#> 7 fquarter4 -0.09295426 0.31180598
#> 8 temp_sc -0.18332144 0.11378328
#> 9 depth_sc -0.31052577 0.14280439
#> 10 density_cod_small_sc -0.11053333 0.22523374
#> 11 density_cod_large_sc 0.16726527 0.19789911
#> 12 density_flounder_sc 0.04175951 0.09689319
#> 13 oxy_sc -0.12249743 0.15996736
#> 14 density_flounder_sc:oxy_sc -0.22168059 0.11837020
sanity(s_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Benthic
s_cod_ben <- sdmTMB(
data = small_cod_stomach,
formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(s_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc + density_cod_large_sc + density_flounder_sc *
#> Formula: oxy_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -5.64 0.42
#> fyear2016 -5.72 0.27
#> fyear2017 -4.97 0.20
#> fyear2018 -5.67 0.21
#> fyear2019 -5.55 0.39
#> fyear2020 -5.90 0.19
#> fquarter4 0.22 0.32
#> temp_sc -0.18 0.12
#> depth_sc -0.36 0.15
#> density_cod_small_sc 0.13 0.23
#> density_cod_large_sc -0.13 0.21
#> density_flounder_sc 0.10 0.10
#> oxy_sc 0.10 0.17
#> density_flounder_sc:oxy_sc -0.29 0.12
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.79
#>
#> Dispersion parameter: 0.14
#> Tweedie p: 1.56
#> ML criterion at convergence: -2943.031
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_ben)
#> term estimate std.error
#> 1 fyear2015 -5.63963437 0.41728201
#> 2 fyear2016 -5.71555073 0.26736376
#> 3 fyear2017 -4.96931045 0.20133219
#> 4 fyear2018 -5.67299371 0.21258020
#> 5 fyear2019 -5.54980300 0.38532858
#> 6 fyear2020 -5.89510007 0.19336637
#> 7 fquarter4 0.21970967 0.32106970
#> 8 temp_sc -0.17986537 0.11810696
#> 9 depth_sc -0.35792181 0.14677782
#> 10 density_cod_small_sc 0.12910755 0.22899007
#> 11 density_cod_large_sc -0.12685493 0.20658268
#> 12 density_flounder_sc 0.09636279 0.09872993
#> 13 oxy_sc 0.09878458 0.16518621
#> 14 density_flounder_sc:oxy_sc -0.29431840 0.12153496
sanity(s_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Saduria
s_cod_sad <- sdmTMB(
data = small_cod_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(s_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc + density_cod_large_sc + density_flounder_sc *
#> Formula: oxy_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -9.08 1.29
#> fyear2016 -8.99 0.84
#> fyear2017 -10.53 0.74
#> fyear2018 -10.22 0.79
#> fyear2019 -10.75 1.33
#> fyear2020 -9.04 0.59
#> fquarter4 0.34 1.06
#> temp_sc 0.14 0.33
#> depth_sc 0.26 0.40
#> density_cod_small_sc -0.15 0.54
#> density_cod_large_sc 0.08 0.55
#> density_flounder_sc -1.34 0.45
#> oxy_sc 0.35 0.46
#> density_flounder_sc:oxy_sc -0.41 0.49
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 2.01
#>
#> Dispersion parameter: 0.19
#> Tweedie p: 1.48
#> ML criterion at convergence: -263.206
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_sad)
#> term estimate std.error
#> 1 fyear2015 -9.07779271 1.2898235
#> 2 fyear2016 -8.98853609 0.8428161
#> 3 fyear2017 -10.52518183 0.7365606
#> 4 fyear2018 -10.21759621 0.7939155
#> 5 fyear2019 -10.75111413 1.3345506
#> 6 fyear2020 -9.04136144 0.5937819
#> 7 fquarter4 0.34186222 1.0616702
#> 8 temp_sc 0.14426222 0.3303548
#> 9 depth_sc 0.26102087 0.4041756
#> 10 density_cod_small_sc -0.15380459 0.5352128
#> 11 density_cod_large_sc 0.08468962 0.5467832
#> 12 density_flounder_sc -1.34276867 0.4522182
#> 13 oxy_sc 0.35260152 0.4586600
#> 14 density_flounder_sc:oxy_sc -0.40751860 0.4911995
sanity(s_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100Large cod
# Total
l_cod_tot <- sdmTMB(
data = large_cod_stomach,
formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(l_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc + density_cod_large_sc + density_flounder_sc *
#> Formula: oxy_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -3.93 0.24
#> fyear2016 -4.49 0.16
#> fyear2017 -4.25 0.13
#> fyear2018 -4.33 0.11
#> fyear2019 -4.44 0.24
#> fyear2020 -4.56 0.13
#> fquarter4 -0.48 0.20
#> temp_sc 0.04 0.08
#> depth_sc 0.05 0.11
#> density_cod_small_sc 0.42 0.22
#> density_cod_large_sc -0.22 0.22
#> density_flounder_sc -0.06 0.06
#> oxy_sc 0.06 0.10
#> density_flounder_sc:oxy_sc -0.03 0.07
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.42
#>
#> Dispersion parameter: 0.60
#> Tweedie p: 1.71
#> ML criterion at convergence: -4812.982
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_tot)
#> term estimate std.error
#> 1 fyear2015 -3.93309844 0.24356798
#> 2 fyear2016 -4.49285009 0.15549512
#> 3 fyear2017 -4.25319947 0.13076704
#> 4 fyear2018 -4.32775670 0.11248077
#> 5 fyear2019 -4.44268256 0.24316977
#> 6 fyear2020 -4.56147977 0.12705887
#> 7 fquarter4 -0.47584719 0.20471977
#> 8 temp_sc 0.03629594 0.07886520
#> 9 depth_sc 0.05126377 0.10663245
#> 10 density_cod_small_sc 0.41917028 0.21975227
#> 11 density_cod_large_sc -0.22222104 0.21800609
#> 12 density_flounder_sc -0.05704968 0.06125595
#> 13 oxy_sc 0.06168315 0.10263081
#> 14 density_flounder_sc:oxy_sc -0.03314841 0.06987393
sanity(l_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Benthic
l_cod_ben <- sdmTMB(
data = large_cod_stomach,
formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(l_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc + density_cod_large_sc + density_flounder_sc *
#> Formula: oxy_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -6.81 0.39
#> fyear2016 -6.65 0.25
#> fyear2017 -6.55 0.21
#> fyear2018 -6.66 0.18
#> fyear2019 -7.78 0.38
#> fyear2020 -6.79 0.20
#> fquarter4 1.07 0.32
#> temp_sc -0.23 0.12
#> depth_sc 0.15 0.17
#> density_cod_small_sc 0.05 0.33
#> density_cod_large_sc 0.29 0.33
#> density_flounder_sc -0.14 0.10
#> oxy_sc 0.34 0.17
#> density_flounder_sc:oxy_sc -0.01 0.12
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.79
#>
#> Dispersion parameter: 0.42
#> Tweedie p: 1.67
#> ML criterion at convergence: -4579.306
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_ben)
#> term estimate std.error
#> 1 fyear2015 -6.810709941 0.39070269
#> 2 fyear2016 -6.650622983 0.25055841
#> 3 fyear2017 -6.548153215 0.20997551
#> 4 fyear2018 -6.659105762 0.18475239
#> 5 fyear2019 -7.778384365 0.38493973
#> 6 fyear2020 -6.786058885 0.20188801
#> 7 fquarter4 1.071428531 0.32284539
#> 8 temp_sc -0.230415833 0.12383263
#> 9 depth_sc 0.149296568 0.17155416
#> 10 density_cod_small_sc 0.047376615 0.33401962
#> 11 density_cod_large_sc 0.285332497 0.33497684
#> 12 density_flounder_sc -0.135768118 0.09772925
#> 13 oxy_sc 0.343186202 0.16712306
#> 14 density_flounder_sc:oxy_sc -0.006315457 0.11754496
sanity(l_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Saduria
l_cod_sad <- sdmTMB(
data = large_cod_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(l_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc + density_cod_large_sc + density_flounder_sc *
#> Formula: oxy_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -8.84 1.26
#> fyear2016 -10.91 0.80
#> fyear2017 -11.11 0.75
#> fyear2018 -11.28 0.67
#> fyear2019 -12.09 1.40
#> fyear2020 -10.54 0.64
#> fquarter4 -0.48 1.04
#> temp_sc 0.25 0.35
#> depth_sc 0.22 0.48
#> density_cod_small_sc 0.03 0.65
#> density_cod_large_sc -0.45 0.71
#> density_flounder_sc -1.08 0.41
#> oxy_sc 0.68 0.47
#> density_flounder_sc:oxy_sc 0.00 0.47
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 2.45
#>
#> Dispersion parameter: 0.12
#> Tweedie p: 1.46
#> ML criterion at convergence: -487.780
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_sad)
#> term estimate std.error
#> 1 fyear2015 -8.843930672 1.2551682
#> 2 fyear2016 -10.910267097 0.8005328
#> 3 fyear2017 -11.107997362 0.7513374
#> 4 fyear2018 -11.278581130 0.6708063
#> 5 fyear2019 -12.089648568 1.4039963
#> 6 fyear2020 -10.540081776 0.6418378
#> 7 fquarter4 -0.477676767 1.0426728
#> 8 temp_sc 0.253550124 0.3549809
#> 9 depth_sc 0.218380287 0.4845376
#> 10 density_cod_small_sc 0.030106804 0.6508573
#> 11 density_cod_large_sc -0.451303087 0.7126573
#> 12 density_flounder_sc -1.083810480 0.4104634
#> 13 oxy_sc 0.675276134 0.4674646
#> 14 density_flounder_sc:oxy_sc 0.004382066 0.4720712
sanity(l_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100Flounder
# Total
fle_tot <- sdmTMB(
data = flounder_stomach,
formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(fle_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc * oxy_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -5.07 0.33
#> fyear2016 -4.96 0.26
#> fyear2017 -4.95 0.17
#> fyear2018 -5.16 0.17
#> fyear2019 -5.01 0.30
#> fyear2020 -4.39 0.14
#> fquarter4 0.47 0.25
#> temp_sc 0.05 0.09
#> depth_sc -0.13 0.12
#> density_cod_small_sc 0.27 0.26
#> oxy_sc 0.21 0.14
#> density_cod_large_sc -0.25 0.26
#> density_flounder_sc -0.08 0.08
#> density_cod_small_sc:oxy_sc -0.13 0.17
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.69
#>
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> ML criterion at convergence: -3635.533
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_tot)
#> term estimate std.error
#> 1 fyear2015 -5.06715201 0.32774948
#> 2 fyear2016 -4.96311363 0.25612380
#> 3 fyear2017 -4.94872102 0.16634274
#> 4 fyear2018 -5.15532320 0.17090355
#> 5 fyear2019 -5.01188780 0.30354649
#> 6 fyear2020 -4.39318789 0.13752021
#> 7 fquarter4 0.47129234 0.25039002
#> 8 temp_sc 0.05425746 0.09198798
#> 9 depth_sc -0.13324303 0.11759711
#> 10 density_cod_small_sc 0.26816835 0.26459501
#> 11 oxy_sc 0.21363168 0.13699214
#> 12 density_cod_large_sc -0.25338400 0.26356414
#> 13 density_flounder_sc -0.08200934 0.07529710
#> 14 density_cod_small_sc:oxy_sc -0.13353163 0.16917662
sanity(fle_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Benthic
fle_ben <- sdmTMB(
data = flounder_stomach,
formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(fle_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc * oxy_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -5.11 0.34
#> fyear2016 -5.00 0.27
#> fyear2017 -4.96 0.17
#> fyear2018 -5.13 0.18
#> fyear2019 -5.30 0.32
#> fyear2020 -4.45 0.14
#> fquarter4 0.47 0.26
#> temp_sc 0.10 0.10
#> depth_sc -0.20 0.12
#> density_cod_small_sc 0.32 0.27
#> oxy_sc 0.23 0.14
#> density_cod_large_sc -0.29 0.27
#> density_flounder_sc -0.06 0.08
#> density_cod_small_sc:oxy_sc -0.09 0.18
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.72
#>
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> ML criterion at convergence: -3574.297
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_ben)
#> term estimate std.error
#> 1 fyear2015 -5.11048157 0.33989780
#> 2 fyear2016 -5.00496148 0.26663702
#> 3 fyear2017 -4.96251227 0.17279852
#> 4 fyear2018 -5.12786504 0.17730557
#> 5 fyear2019 -5.29751455 0.31926186
#> 6 fyear2020 -4.44609046 0.14312220
#> 7 fquarter4 0.46914758 0.26092314
#> 8 temp_sc 0.09771251 0.09572682
#> 9 depth_sc -0.19742836 0.12302319
#> 10 density_cod_small_sc 0.31634925 0.27402938
#> 11 oxy_sc 0.22536646 0.14331473
#> 12 density_cod_large_sc -0.28796091 0.27286263
#> 13 density_flounder_sc -0.05528865 0.07821425
#> 14 density_cod_small_sc:oxy_sc -0.08752974 0.17644492
sanity(fle_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Saduria
fle_sad <- sdmTMB(
data = flounder_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(fle_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: density_cod_small_sc * oxy_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: tweedie(link = 'log')
#>
#> coef.est coef.se
#> fyear2015 -6.33 1.21
#> fyear2016 -5.62 0.91
#> fyear2017 -8.61 0.61
#> fyear2018 -10.12 0.69
#> fyear2019 -8.55 1.18
#> fyear2020 -8.48 0.56
#> fquarter4 -1.26 0.91
#> temp_sc -0.40 0.31
#> depth_sc 0.03 0.37
#> density_cod_small_sc -0.98 0.63
#> oxy_sc 0.26 0.43
#> density_cod_large_sc 0.51 0.62
#> density_flounder_sc -1.00 0.37
#> density_cod_small_sc:oxy_sc -0.35 0.56
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 2.63
#>
#> Dispersion parameter: 0.21
#> Tweedie p: 1.49
#> ML criterion at convergence: -1086.186
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_sad)
#> term estimate std.error
#> 1 fyear2015 -6.32974501 1.2138097
#> 2 fyear2016 -5.62234282 0.9086070
#> 3 fyear2017 -8.60780551 0.6072356
#> 4 fyear2018 -10.12096544 0.6889580
#> 5 fyear2019 -8.55498870 1.1832802
#> 6 fyear2020 -8.48135360 0.5636020
#> 7 fquarter4 -1.25920836 0.9118550
#> 8 temp_sc -0.40038843 0.3143738
#> 9 depth_sc 0.03433966 0.3686257
#> 10 density_cod_small_sc -0.98458983 0.6321219
#> 11 oxy_sc 0.26364504 0.4285877
#> 12 density_cod_large_sc 0.50812859 0.6226676
#> 13 density_flounder_sc -1.00034650 0.3727599
#> 14 density_cod_small_sc:oxy_sc -0.35448644 0.5617143
sanity(fle_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100# mcmc_res_l_cod_sad_v2 <- residuals(l_cod_sad_v2,
# type = "mle-mcmc",
# mcmc_iter = 201,
# mcmc_warmup = 200)
#
# qqnorm(mcmc_res_l_cod_sad_v2, ylim = c(-5, 5), xlim = c(-5, 5)); qqline(mcmc_res_l_cod_sad_v2)Plot coefficients
#### Total prey
# small cod total prey
small_cod_tot_pars <- bind_rows(tidy(s_cod_tot, effects = "ran_pars", conf.int = TRUE),
tidy(s_cod_tot, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Small cod",
response = "Total prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# large cod total prey
large_cod_tot_pars <- bind_rows(tidy(l_cod_tot, effects = "ran_pars", conf.int = TRUE),
tidy(l_cod_tot, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Large cod",
response = "Total prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# flounder total prey
flounder_tot_pars <- bind_rows(tidy(fle_tot, effects = "ran_pars", conf.int = TRUE),
tidy(fle_tot, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Flounder",
response = "Total prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
#### Benthic prey
# small cod benthic prey
small_cod_ben_pars <- bind_rows(tidy(s_cod_ben, effects = "ran_pars", conf.int = TRUE),
tidy(s_cod_ben, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Small cod",
response = "Benthic prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# large cod benthic prey
large_cod_ben_pars <- bind_rows(tidy(l_cod_ben, effects = "ran_pars", conf.int = TRUE),
tidy(l_cod_ben, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Large cod",
response = "Benthic prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# flounder benthic prey
flounder_ben_pars <- bind_rows(tidy(fle_ben, effects = "ran_pars", conf.int = TRUE),
tidy(fle_ben, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Flounder",
response = "Benthic prey FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
#### Saduria
# small cod saduria
small_cod_sad_pars <- bind_rows(tidy(s_cod_sad, effects = "ran_pars", conf.int = TRUE),
tidy(s_cod_sad, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Small cod",
response = "Saduria FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# large cod saduria
large_cod_sad_pars <- bind_rows(tidy(l_cod_sad, effects = "ran_pars", conf.int = TRUE),
tidy(l_cod_sad, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Large cod",
response = "Saduria FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# flounder saduria
flounder_sad_pars <- bind_rows(tidy(fle_sad, effects = "ran_pars", conf.int = TRUE),
tidy(fle_sad, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Flounder",
response = "Saduria FR")
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
coefs <- bind_rows(small_cod_tot_pars,
large_cod_tot_pars,
flounder_tot_pars,
small_cod_ben_pars,
large_cod_ben_pars,
flounder_ben_pars,
small_cod_sad_pars,
large_cod_sad_pars,
flounder_sad_pars)
coefs <- coefs %>%
mutate(term2 = term,
term2 = ifelse(term == "temp_sc", "Temperature", term2),
term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
term2 = ifelse(term == "depth_sc", "Depth", term2),
term2 = ifelse(term == "density_cod_large_sc", "Cod density (> 29)", term2),
term2 = ifelse(term == "density_cod_small_sc", "Cod density (< 29)", term2),
term2 = ifelse(term == "density_cod_small_sc:oxy_sc", "Cod density (> 29) x Oxygen", term2),
term2 = ifelse(term == "density_flounder_sc:oxy_sc", "Flounder density x Oxygen", term2),
term2 = ifelse(term == "density_cod_small_sc:oxy_sc", "Small cod density x Oxygen", term2),
term2 = ifelse(term == "density_flounder_sc", "Flounder density", term2))
#> mutate: new variable 'term2' (character) with 18 unique values and 0% NA
# Fixed continuous effects:
coefs %>%
filter(!term %in% c("phi", "sigma_G", "tweedie_p")) %>%
filter(!grepl('year', term)) %>%
filter(!grepl('quarter', term)) %>%
ggplot(aes(term2, estimate, color = species)) +
geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
geom_point(position = position_dodge(width = 0.8)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
position = position_dodge(width = 0.8), alpha = 0.5) +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
facet_wrap(~response) +
labs(x = "Predictor", y = "Standardized coefficient") +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
aspect.ratio = 1.5) +
NULL
#> filter: removed 27 rows (18%), 126 rows remaining
#> filter: removed 54 rows (43%), 72 rows remaining
#> filter: removed 9 rows (12%), 63 rows remaining
ggsave("figures/fr_coefs.pdf", width = 17, height = 17, units = "cm")
coefs %>%
filter(term %in% c("fyear2015", "fyear2016", "fyear2017", "fyear2018", "fyear2019", "fyear2020", "quarter")) %>%
mutate(term2 = ifelse(term == "fyear2015", "2015", term2),
term2 = ifelse(term == "fyear2016", "2016", term2),
term2 = ifelse(term == "fyear2017", "2017", term2),
term2 = ifelse(term == "fyear2018", "2018", term2),
term2 = ifelse(term == "fyear2019", "2019", term2),
term2 = ifelse(term == "fyear2020", "2020", term2),
term2 = ifelse(term == "quarter", "Quarter", term2)) %>%
filter(!term2 == "Quarter") %>%
ggplot(aes(term2, estimate, color = species)) +
geom_point(position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
facet_wrap(~response, scales = "free") +
labs(x = "Predictor", y = "Standardized coefficient") +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
aspect.ratio = 1.5) +
NULL
#> filter: removed 99 rows (65%), 54 rows remaining
#> mutate: changed 54 values (100%) of 'term2' (0 new NA)
#> filter: no rows removed
ggsave("figures/supp/fr_coefs_yr.pdf", width = 17, height = 17, units = "cm")# fle_sad <- sdmTMB(
# data = flounder_stomach,
# formula = saduria_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = flounder_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
s_cod_sad2 <- sdmTMB(
data = small_cod_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
l_cod_sad2 <- sdmTMB(
data = large_cod_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
fle_sad2 <- sdmTMB(
data = flounder_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc,
family = tweedie(link = "log"),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 7)),
c(rep(NA, 7), rep(1, 7)))),
control = sdmTMBcontrol(newton_loops = 1)
)
quantile(flounder_stomach$density_flounder, prob = c(0.1, 0.9))
#> 10% 90%
#> 231.5209 4716.5631
nd_fle <- data.frame(density_flounder = seq(0, 5000, length.out = 100))
nd_fle <- nd_fle %>%
mutate(density_flounder_sc = (density_flounder - mean(flounder_stomach$density_flounder)) / sd(flounder_stomach$density_flounder),
year = 2017L,
fyear = as.factor(2017),
fquarter = as.factor(1),
temp_sc = 0,
depth_sc = 0,
oxy_sc = 0,
density_cod_small_sc = 0,
density_cod_large_sc = 0
)
#> mutate: new variable 'density_flounder_sc' (double) with 100 unique values and 0% NA
#> new variable 'year' (integer) with one unique value and 0% NA
#> new variable 'fyear' (factor) with one unique value and 0% NA
#> new variable 'fquarter' (factor) with one unique value and 0% NA
#> new variable 'temp_sc' (double) with one unique value and 0% NA
#> new variable 'depth_sc' (double) with one unique value and 0% NA
#> new variable 'oxy_sc' (double) with one unique value and 0% NA
#> new variable 'density_cod_small_sc' (double) with one unique value and 0% NA
#> new variable 'density_cod_large_sc' (double) with one unique value and 0% NA
# Predict from full model
sc_margin_sad <- predict(s_cod_sad2, newdata = nd_fle, se_fit = TRUE, re_form = NA, re_form_iid = NA) %>%
mutate(species2 = "Small cod")
#> mutate: new variable 'species2' (character) with one unique value and 0% NA
lc_margin_sad <- predict(l_cod_sad2, newdata = nd_fle, se_fit = TRUE, re_form = NA, re_form_iid = NA) %>%
mutate(species2 = "Large cod")
#> mutate: new variable 'species2' (character) with one unique value and 0% NA
f_margin_sad <- predict(fle_sad2, newdata = nd_fle, se_fit = TRUE, re_form = NA, re_form_iid = NA) %>%
mutate(species2 = "Flounder")
#> mutate: new variable 'species2' (character) with one unique value and 0% NA
cond_effects <- bind_rows(sc_margin_sad, lc_margin_sad, f_margin_sad)
ggplot(cond_effects, aes(density_flounder, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se),
color = species2,
fill = species2,
linetype = species2)) +
geom_ribbon(alpha = 0.2, color = NA) +
theme(aspect.ratio = 1) +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_fill_brewer(palette = "Paired", name = "Species", direction = -1) +
scale_linetype_manual(values = c(3, 2, 1)) +
guides(linetype = "none",
color = guide_legend(override.aes = list(linetype = c(3, 2, 1)))) +
geom_line(size = 1.3) +
labs(x = "Flounder biomass density",
y = "Saduria feeding ratio") +
theme(legend.position = c(0.9, 0.9))
ggsave("figures/flounder_saduria_conditional.pdf", width = 17, height = 17, units = "cm")